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Using  the  multiphase  free-energy  lattice  Boltzmann  method  (LBM),  the  formation  of  a  water  droplet 
emerging  through  a  micro-pore  on  the  hydrophobic  gas  diffusion  layer  (GDL)  surface  in  a  proton  exchange 
membrane  fuel  cell  (PEMFC)  and  its  subsequent  movement  on  the  GDL  surface  under  the  action  of  gas 
shear  are  simulated.  The  dynamic  behavior  of  the  water  droplet  emergence,  growth,  detachment  and 
movement  in  the  gas  flow  channel  is  presented.  The  size  of  the  detached  droplet  and  the  time  of  the 
droplet  removing  out  of  the  channel  under  the  influence  of  gas  flow  velocity  and  GDL  surface  wettability 
are  investigated.  The  results  show  that  water  droplet  removal  is  facilitated  by  a  high  gas  flow  velocity 
on  a  more  hydrophobic  GDL  surface.  A  highly  hydrophobic  surface  is  shown  to  be  capable  of  lifting  the 
water  droplet  from  the  GDL  surface,  resulting  in  more  GDL  surface  available  for  gas  reactant  transport. 
Furthermore,  an  analytical  model  based  on  force  balance  is  presented  to  predict  the  droplet  detachment 
size,  and  the  predicted  results  are  in  good  agreement  with  the  simulation  results.  It  is  shown  that  the 
LBM  approach  is  an  effective  tool  to  investigate  water  transport  phenomena  in  the  gas  flow  channel  of 
PEMFCs  with  surface  wettability  taken  into  consideration. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Proton  exchange  membrane  fuel  cells  (PEMFCs),  as  an  attractive 
power  source  for  automotive  and  portable  applications,  have  been 
given  much  attention  in  the  last  decade.  In  a  PEMFC  cathode,  oxy¬ 
gen  is  supplied  to  the  gas  flow  channel  and  transfers  to  the  catalyst 
layer  where  electrochemical  reactions  occur  and  water  is  produced 
as  shown  in  Fig.  1(a).  The  transport  of  the  produced  water  has  a 
significant  influence  on  the  performance  of  PEMFCs  [1  ].  In  order  to 
maintain  high  proton  conductivity,  it  is  necessary  to  keep  sufficient 
water  content  in  the  membrane.  On  the  other  hand,  too  much  water 
in  the  cathode  can  cause  flooding  in  the  porous  gas  diffusion  layer 
(GDL)  or  clogging  of  the  gas  flow  channel  by  the  accumulated  liq¬ 
uid  water,  thus  resulting  in  the  inhibition  of  the  gaseous  reaction 
transport  to  the  reaction  sites  [2].  Therefore,  water  management  is 
considered  to  be  a  critical  issue  on  the  performance  of  PEMFCs. 

Previous  studies  on  the  water  transport  in  the  PEMFCs  have 
mainly  focused  on  the  water  distribution  in  the  porous  GDL,  which 
has  been  assumed  to  be  a  homogeneous  medium  [3-6].  Although 
the  clogging  of  the  gas  channel  could  cause  severe  increase  of  the 
gas  flow  pressure  loss  thus  leading  to  more  serious  degradation  of 
fuel  cell  performance  than  GDL  flooding,  water  transport  in  the  gas 
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flow  channel  has  not  been  given  adequate  attention  until  recently. 
In  almost  all  of  the  numerical  studies,  water  in  the  gas  channel  is 
simply  assumed  to  be  vapor  flow  or  mist  flow  [7].  Visualization 
studies,  however,  did  not  support  such  an  assumption.  In  practice, 
water  is  often  accumulated  to  form  individual  droplets  through 
pores  emerging  on  the  GDL  surface  [8-11],  which  is  schematically 
shown  in  Fig.  1(a).  Yang  et  al.  [8]  performed  a  detail  optical  visu¬ 
alization  experiment  on  water  droplet  dynamics  at  high  current 
densities  and  operating  temperatures  representative  of  automotive 
PEMFCs,  and  confirmed  the  appearance  of  water  droplets  in  the  gas 
flow  channel.  Using  high-resolution  images,  the  mechanics  of  liquid 
water  emergence  on  the  hydrophobic  GDL  surface,  droplet  growth 
and  departure,  as  well  as  the  two-phase  flow  in  the  gas  flow  channel, 
is  characterized.  Recently,  neutrons  imaging  technique  [9,10]  and 
gas  chromatography  measurements  [11],  as  well  as  pressure  drop 
measurement  [12,13]  which  is  related  to  the  water  content  in  the 
gas  flow  channel,  have  also  been  used  to  probe  water  droplet  forma¬ 
tion  and  distribution  in  the  gas  flow  channel  in  operating  PEMFCs, 
and  similar  visualization  results  have  been  shown.  These  studies 
have  provided  some  basic  understanding  of  water  transport  in  the 
gas  flow  channel  but  they  were  all  limited  to  qualitative  investi¬ 
gations.  It  has  been  known  that  the  size  and  distribution  of  the 
water  droplet,  as  well  as  the  geometry  of  the  gas  channel  and  the 
wettability  property  of  the  GDL  have  remarkable  effects  on  water 
removal.  Thus,  in  order  to  obtain  a  thorough  understanding  of  this 
issue,  quantitative  results  are  needed. 
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Nomenclature 

A,  B  free-energy  parameter 

Bo  bond  number 

b  width  of  the  channel  above  the  droplet 

Ca  capillary  number 

c  lattice  speed 

cs  sound  speed 

D  droplet  diameter 

d  GDL  pore  equivalent  diameter 

e  lattice  velocity  vector 

F  body  force  vector 

Fp  pressure  force 

Fshear  shear  force 

F surface  surface  force 

/,  g  distribution  function 

H  channel  height 

h  droplet  height 

1  unit  tensor 

L  contact  line  length 

l0  length  scale 

M  mobility 

P  pressure  tensor 

p  pressure 

Qd  water  flow  rate 

R  droplet  radius 

t0  time  scale 

tm  droplet  moving  time 

tf  droplet  formation  time 

U  velocity  magnitude 

u  velocity  vector 

Vd  detached  droplet  volume 

w  weight  coefficient 

x  position  vector 

Greek  symbols 

r  mobility  coefficient 

8t  time  step 

8x  lattice  distance 

r\  dynamic  viscosity 

0  wetting  potential 

6  contact  angle 

k  surface  tension  parameter 

fi  chemical  potential 

§  phase  interface  width 

p  density 

a  surface  tension 

tf,  Tg  relaxation  parameter 
v  kinematic  viscosity 

order  parameter 

\j/  bulk  free-energy  density 

Superscripts 

a  air 

ch  chemical 

d  droplet 

eq  equilibrium 

/  fluid 

by  hydrostatic 

in  inlet 

k  dimension 

l  lattice  unit 

p  physical  unit 

s  solid 

th  thermodynamic 

a  ath  direction  of  the  lattice  velocity 


More  recently,  Zhang  et  al.  [14]  analyzed  the  forces  acting  on 
a  static  water  droplet  on  the  GDL  surface  based  on  their  optical 
observations  of  an  operating  fuel  cell.  The  predicted  detachment 
diameter  of  water  droplets  emerging  on  the  GDL  was  compared 
with  experimental  results.  Chen  et  al.  [15]  proposed  a  simplified 
model  to  predict  the  onset  of  water  droplet  instability  on  a  GDL 
under  the  gas  flow  shear  based  on  a  macroscopic  force  balance 
analysis.  The  gas  flow  velocities,  channel  heights  and  contact  angle 
hysteresis  (the  difference  between  the  advancing  and  receding  con¬ 
tact  angles)  were  taken  into  consideration  to  characterize  droplet 
instability.  The  droplet  instability  windows  predicted  by  this  model 
were  compared  with  visualization  experiments  and  2D  finite  ele¬ 
ment  simulations.  Following  Chen’s  work  [15],  Kumbur  et  al.  [16] 
proposed  a  similar  analytical  model  to  study  effects  of  channel 
geometry  and  GDL  surface  properties  such  as  PTFE  contents  on 
water  droplet  instability.  The  contact  angle  hysteresis  was  directly 
introduced  into  their  model  as  an  empirical  correlation  of  PTFE  con¬ 
tent  from  their  experiment  results  rather  than  assuming  a  constant 
value.  However,  these  works  did  not  provide  sufficient  quantita¬ 
tive  understanding  of  water  transport  on  the  GDL  surface  as  they 
were  focused  on  a  static  water  droplet,  which  is  different  from  the 
observed  dynamic  characteristic  of  water  droplet  growth,  departure 
and  subsequent  movement.  At  the  same  time,  Kimball  et  al.  [17,18] 
and  Benziger  et  al.  [19]  presented  a  different  mechanism  of  water 
droplet  removal  in  gas  flow  channel,  where  the  emerging  droplets 
holding  on  the  GDL  surface  is  attributed  to  the  surface  energy 
of  the  water  and  the  GDL  pore  but  not  the  interaction  between 
the  droplet  and  GDL  surface.  Additionally,  water  droplet  dynamic 
behavior  on  solid  surfaces  is  a  fundamental  fluid  mechanics  prob¬ 
lem  with  important  applications  not  only  in  fuel  cells  but  also  in 
many  other  systems  such  as  chemical  process  technologies  [15]  as 
well.  In  spite  of  this,  the  problem  has  not  been  thorough  inves¬ 
tigated  theoretically  because  of  the  difficulty  associated  with  the 
contact  line  dynamics,  which  is  still  a  pending  physical  problem  to 
be  dealt  with  to-date  [20].  Several  CFD  studies  based  on  volume- 
of-fluid  (VOF)  method  have  also  been  performed  to  simulate  the 
water  droplet  dynamic  behavior  and  distribution  in  a  serpentine 
channel  with  an  initially  prescribed  water  volume  in  the  chan¬ 
nel  [21,22],  but  the  local  wettability  property  of  the  GDL  surface 
and  the  water  emergence  process  were  not  explicitly  considered. 
Although  the  effects  of  gas  flow  velocity  and  GDL  wettability  on 
water  droplet  deformation  [23,24]  and  formation  [25,26]  have  also 
been  studied  most  recently  by  VOF  model,  the  VOF  method  is  defi¬ 
cient  to  model  dynamic  contact  angle  change,  which  is  a  critical 
parameter  in  the  droplet  dynamic  process  in  the  present  problem. 
In  practice,  the  interface  topology  in  the  VOF  method  is  artificially 
reconstructed  by  calculating  macroscopic  curvature,  which  has  lit¬ 
tle  physical  basis  [27,28].  Usually,  in  order  to  obtain  the  dynamic 
change  of  the  contact  angle,  a  complicated  numerical  scheme  must 
be  used  to  track  interface  change  continuously  in  the  VOF  method 
[20].  Strictly  speaking,  the  interface  between  different  phase  and 
the  contact  line  dynamics  on  the  solid  surface  were  based  on  the 
mesoscopic  scale  [29].  Thus,  mesoscopic  level  studies  are  expected 
to  investigate  accurately  the  effect  of  surface  wettability  of  the  GDL 
on  water  droplet  dynamic  behaviors  in  the  gas  channel. 

On  the  other  hand,  the  lattice  Boltzmann  method  (LBM),  based 
on  the  mesoscale  kinetic  theory  of  fluids,  has  been  developed  and 
applied  to  simulate  fluid  flow  and  various  transport  problems  in 
recent  years.  In  particular,  due  to  the  kinetic  nature  and  no  need 
to  track  the  phase  interface  explicitly,  the  LBM  has  been  found  to 
be  a  very  effective  numerical  tool  to  simulate  multiphase,  multi- 
component  flow  problems  [30-32].  More  recently,  some  works  on 
the  application  of  the  LBM  to  study  the  transport  phenomena  in 
PEMFCs  have  been  carried  out  [33-37].  For  example,  Fei  and  Hong 
[33]  simulated  the  two-phase  flow  of  methanol/C02  in  a  micro- 
channel  as  a  flow  passage  in  the  GDL  of  a  micro-direct  methanol  fuel 
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Fig.  1.  (a)  Schematic  of  the  transport  process  in  a  PEMFC  cathode,  (b)  Schematic  of  the  force  analysis  of  an  emerging  droplet. 


cell,  and  studied  CO2  bubbly  flow  phenomena  in  the  micro-channel 
under  different  operation  conditions.  Joshi  et  al.  [34]  used  LBM  to 
model  the  mass  transport  of  H2  and  the  produced  H20  (vapor)  in 
the  presence  of  N2  in  the  2D  porous  anode  structure  of  a  solid  oxide 
fuel  cell  (SOFC).  Suzue  et  al.  [35]  also  performed  a  LBM  simulation 
on  multicomponent  flow  in  a  SOFC.  Park  and  Li  [36]  presented  a 
2D  two-phase  LBM  simulation  of  water  droplet  flowing  through 
a  fibrous  structure  of  carbon  paper  GDL.  Furthermore,  Niu  et  al. 
[37]  applied  the  LBM  model  to  simulate  two-phase  transport  in 
the  carbon  paper  GDL  of  a  PEMFC,  and  investigated  the  relative 
permeability  under  different  conditions.  Flowever,  the  application 
of  the  LBM  to  study  water  droplet  dynamics  on  the  hydrophobic 
GDL  surface  adjacent  to  a  gas  channel  has  not  yet  been  performed 
to-date. 

In  this  work,  the  effects  of  gas  flow  velocity  and  GDL  wettability 
on  water  droplet  dynamic  behavior  including  emergence,  deforma¬ 
tion,  departure  and  subsequent  movement  in  the  gas  flow  channel 
with  a  hydrophobic  GDL  surface  are  simulated  using  the  multiphase 
LBM  approach  based  on  the  free-energy  theory  [38].  Fig.  1(a)  shows 
a  schematic  of  the  actual  water  emergence  process  with  several 
droplet  formation  spots  on  a  GDL  surface  as  visualized  in  previ¬ 
ous  experiments,  and  Fig.  1(b)  shows  its  simplified  version  of  only 
one  water  emergence  pore,  which  will  be  numerically  simulated 
in  this  paper.  In  addition,  based  on  the  forces  acting  on  the  droplet 
shown  in  Fig.  1  (b),  a  simple  analytical  model  is  developed  to  predict 
the  droplet  detachment  size,  and  comparison  is  made  between  the 
model  predictions  and  the  simulation  results. 


2.  Description  of  the  simulation  model 

At  the  present  time,  the  intermolecular  potential  model  and 
free-energy  model  [39]  are  two  major  approaches  used  in  the  LBM 
to  describe  multiphase  systems.  The  former  model  introduces  the 
nearest-neighboring  interaction  between  fluid  particles  to  describe 
the  intermolecular  potential,  and  the  phase  separation  occurs  with 
a  proper  potential  function  be  chosen.  Flowever,  as  pointed  out  by 


He  and  Doolen  [39],  the  surface  tension  in  this  model  is  actually 
a  numerical  artifact  and  its  value  is  controlled  by  a  force  param¬ 
eter,  which  cannot  be  prescribed  a  priori.  In  the  latter  model,  the 
free-energy  function  is  incorporated  into  the  LBM  to  describe  mul¬ 
tiphase  or  multicomponent  system,  which  has  a  firm  theoretical 
basis.  The  surface  tension  as  well  as  the  contact  angle  when  solid 
is  present  can  be  directly  derived  from  the  free-energy  function 
based  on  Cahn-Hilliard  theory  [38].  Due  to  its  more  physical  basis 
than  intermolecular  potential  model,  the  free-energy  LBM  model 
has  been  widely  used  to  simulate  droplet  formation,  deformation 
and  breakup  [32,40].  Thus,  this  method  is  adopted  in  the  present 
study. 

2.1.  LBM  model 


In  the  free-energy  LBM  model  for  two-phase  system  such  as  gas 
and  liquid  which  have  density  of  pg  and  p/,  respectively,  two  dis¬ 
tribution  functions  /a(x,  t)  and  ga(x,  t)  are  used  to  model  density 
p  =  pg  +  pi ,  velocity  u  and  order  parameter  cp  which  denotes  the 
different  phases,  respectively.  The  evolutions  of  the  distribution 
functions  are  described  by  the  following  equations: 


fair +  ea8t,t  +  8t)~  fa(r,  t ) 

=  ^r\faq(r  +  ea8t,  t  +  8t)-fa{r,  t)]  +  Fa  (la) 

rf 

ga{r+ea8t,  t+8t)-ga{r,  t)=^[geaq{r+ea8t,  t+8t)-ga{r,  t)]  (lb) 

Tg 

where  z y  and  rg  are  two  independent  relaxation  parameters;  ea  is 
lattice  velocity  in  the  ath  direction  and  the  corresponding  lattice 
velocity  directions  of  three-dimension  19-velocity  (D3Q19)  model 
are  shown  in  Fig.  2 .f^q  and geq  are  the  equilibrium  distributions  of fa 
and  ga ,  which  are  given  as  the  functions  of  local  density  p,  velocity 
u  and  order  parameter  cp,  respectively: 


fo^  —  Aa  +  Wap 


e<*  «  ( ea-u)2 

c?  2c? 


u  u 

2 W 


(2a) 
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Fig.  2.  Lattice  velocity  directions  of  the  D3Q19  LBM. 


ga  =Ha+  Wa(p 


u  ( ea-uf 
2c? 


u  u 

2c[ 


(2b) 


where  Aa  and  Ha  are  the  indeterminate  coefficients  and  cs  denotes 
the  speed  of  sound  which  is  given  by  c/V 3  with  c  =  8x/8t  being  the 
lattice  speed  and  8x  being  the  lattice  distance. 

For  the  D3Q19  model,  the  lattice  velocity  vector  ea  and  the 
weight  coefficients  wa  which  depend  on  the  lattice  velocity  direc¬ 
tions  are  given  as  follows  [30]: 

r  (o,o,o), 

ea=  1  (±1,0,0 )c,  (0,  ±l,0)c, 

[(±l,±l,0)c,  (±1,0,  ±l)c, 

{1/3,  a  =  0;  "j 

1/18,  a  =  l,2,...,6;  > 

1/36,  a  =  7,  8,...,  18.  J 

Fa  in  Eq.  ( 1  a)  represents  the  body  force  components  in  lattice  space, 
which  is  given  by  [41  ], 


a  -  0; 

(0,  0  ±  l)c,  a  =  1,2,...,  6; 
(0,  ±1,  ±l)c,  a  =  7,8,  ...,18. 


(3) 


Fa  —  Wa 


u  ea  u 

I  a  “ol 


F 


(4) 


with  F  the  body  force  such  as  gravity. 

The  macroscopic  local  density  p(x,  t),  momentum  j(x,  t)  and 
order  parameter  (p(x,  t)  are  related  to  the  distribution  functions /a(x, 
t)  and  ga  (x,  t)  by 


18 


18 


P=J>  j  =  pu  =  y faea  +  jF;  <p  =  ^gq 


18 


(5) 


a=0 


a=0 


a=0 


Since  these  macroscopic  quantities  must  be  locally  conserved 
in  any  collision  process,  Eq.  (5)  should  also  be  satisfied  by  equi¬ 
librium  distribution  functions/^9  and  g^9.  Moreover,  the  resulting 
continuum  equations  from  the  LBM  model  must  exactly  describe 
the  dynamics  of  a  two-phase  system  [38].  Therefore,  additional 
higher  moments  of f^q  and  g^9  must  satisfy, 


'fffaqeael  =  Pth  +  puuT;  ffglqea  =  <pu\ 


a=0 


a=0 


Tjffeael  =  r ul  +  <puuT 


a=0 


(6) 


where  Pth  is  the  complete  thermodynamic  pressure  tensor,  r  is 
the  coefficient  which  controls  the  phase  interface  diffusion  and  is 
related  to  the  mobility  M  of  the  fluid  as  follows: 


M  =  r(rg  -  0.5)8t 


(7) 


The  pressure  tensor  Pth  and  the  surface  tension  in  two-phase 
system,  as  well  as  the  wetting  boundary  condition  with  the  pres¬ 
ence  of  solid  wall  can  be  determined  from  the  free-energy  basis  as 
below. 


2.2.  Free-energy  model 


For  a  binary  liquid  system  with  the  presence  of  solid  boundary, 
the  free-energy  function  F  can  be  described  as  a  function  of  the 
order  parameter  (p  as  follows  [42]: 

F  =  J  dV  (W)  +  V^)2)  +  jdS  f{(ps)  (8) 

where  the  first  part  on  the  right-hand  side  is  the  Landau  free-energy 
function,  describing  equilibrium  properties  of  the  fluid  system  in 
the  bulk  control  volume  V.  \ fr{(p)  is  the  bulk  free-energy  density  for 
the  homogeneous  system,  and  cps  is  the  wall  order  parameter  which 
is  used  to  denote  the  wall  wettability.  The  square  of  the  gradient 
term  denotes  the  free  energy  excess  in  the  interfacial  region,  which 
is  defined  as  the  surface  energy  between  different  fluids  with  /c 
relating  to  the  surface  tension.  The  free-energy  density  \ /s((p)  can  be 
chosen  as  the  form  of  a  double-well  potential  [40]: 

(9) 

where  A  and  B  are  two  parameters,  which  are  chosen  by  positive 
values  for  phase  separation  in  present  study.  From  the  Landau  free- 
energy  function,  the  chemical  potential  fi  is  given  by  [40]: 

li  =  -A(p  +  Bcp3  -  K?V2(p  (10) 

Following  the  Gibbs-Duhem  equation,  it  gives: 

S/-Pch=(pVp  (11) 

with  the  chemical  part  of  the  thermodynamic  pressure  tensor  Pdl 
given  by 

pch  =p0i  +  K{v<p)(y<p)1  (12) 

where  1  is  the  second-order  unit  tensor,  p0  is  the  scale  part  of  the 
pressure  tensor  which  can  be  expressed  by  [32,38]: 

Po=(pu  -  Vd»-|(V<p)2=  -  ^<P2+^-<P2  -  K<pV2y>-|(Vy>)2  (13) 

The  complete  thermodynamic  pressure  tensor  Pth  is  the  sum  of 
Pch  and  the  hydrostatic  pressure  phy  as  follows: 

Pth=Pch+phyl  (14) 

where  phy  is  given  by  phy  =  c]p  with  cs  being  the  sound  speed  in 
the  Lattice  Boltzmann  scheme.  From  Eqs.  (12)-(14),  the  final  bulk 
pressure  of  the  binary  liquid  system  can  be  obtained  as 

p  =  c2p-^<p2  +  ^<p4-/cy>V2y>-|(V<p)2  (15) 

The  minimization  of  the  derivation  of  the  free  energy  in  bulk 
with  respect  to  order  parameter  cp  results  in  [40]: 

Li  = -A(p  +  B(p3  -  kV2 (p  =  0  (16) 


Thus,  the  equilibrium  order  parameter  <p0  in  the  bulk  of  the 
two  different  fluids  follows  (po  =  i^/B)1/2  =  ±1  when  restricting  the 
value  A  =  B. 

In  addition,  for  a  plane  phase  interface  under  equilibrium  con¬ 
ditions,  the  order  parameter  profile  across  the  phase  interface  can 
be  represented  as, 


(p  =  <p0tanh 


(17) 
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where  £  is  the  interface  width  which  is  given  by  §  =  (2 /c/A)1/2  [32]. 
The  interface  surface  tension  is  evaluated  according  to  thermody¬ 
namic  theory  by  [38], 


From  Eqs.  (17)  and  (18),  the  surface  tension  is  obtained  by, 


4  K  n 


3? 


(19) 


When  specifying  the  surface  tension  o  and  interface  width  §, 
the  chemical  potential  /x  and  thermodynamic  pressure  tensor  Pfh, 
which  are  used  in  the  free-energy  LBM  model  above,  can  be  derived. 
The  interface  width  §  is  a  free  parameter  in  this  model.  In  order  to 
keep  a  sharp  phase  interface,  the  interface  width  should  be  cho¬ 
sen  with  a  small  value.  However,  too  small  §  will  lead  to  numerical 
instability  and  inaccuracy.  In  this  study,  we  specified  its  value  as  two 
lattice  distance,  which  has  been  proved  to  give  stable  and  accuracy 
numerical  results  [40]. 

Applying  Eq.  (14),  the  coefficients  Aa  and  Ha  in  the  equilibrium 
distribution  functions  given  by  Eqs.  (2a)  and  (2b)  can  be  determined 
from  Eqs.  (5)  and  (6)  as  follows: 


'  p-2(P&  +  P!»+Pg)l3c2  a  =  0; 

Aa  =  i  eTaPtheal4c4  -  (P*  +  P$  +  P®  )/36c2  a  =  1,2,... 
,  e£Pthea/8c4  -  (P%  +  P$  +  Pg )/18c2  a  =  7,  8, . . . 

{(p  -  2/7 x/c2  a  =  0; 

T/X/6C2 


6; 
18. 


T/x/12c2 


a  =  1,2,  ...,6; 

a  =  7,  8,  ...,18.  J 


> 


(20) 


with  A  being  a  constant  depending  on  the  contact  angle  as  will 
be  discussed  later.  As  a  result,  minimizing  the  free-energy  func¬ 
tion  F  variations  at  equilibrium  condition  results  in  the  following 
expression  near  the  wall: 


X  =  K 


dip 

dn 


(24) 


where  h  denotes  the  direction  normal  to  the  wall  surface  and  point¬ 
ing  into  the  liquid.  The  parameter  X  can  further  be  obtained  by 
[42]: 


with  0  being  the  wetting  potential,  which  is  related  to  the  contact 
angle  as  follows: 


cos  (0)  = 


(i + @)3/2  -  (i  -  @)3/2 

2 


(26) 


From  Eq.  (26),  the  wetting  potential  can  also  be  obtained  by  equi¬ 
librium  contact  angle  as  [43]: 


0  = 


=2sign(f-p) 


COS'-  I  1  -  cos 


1/2 


(27) 


where  /3  =  arccos(sin2  6)  and  sign(7r/2  -  6)  is  the  sign  function  with 
0  being  the  prescribed  contact  angle.  In  practice,  when  an  equilib¬ 
rium  contact  angle  is  prescribed,  the  wetting  potential  0  and  the 
parameter  X  are  derived  from  Eqs.  (27)  and  (25),  respectively.  As 
a  result,  the  order  parameter  derivative  dip/dn  at  the  wall  is  deter¬ 
mined  from  Eq.  (24),  which  can  be  regarded  as  the  wall  boundary 
condition  for  order  parameter  ip. 


3.  Boundary  conditions 


Employing  the  Chapman-Enskog  multiscale  analysis  for  D3Q19 
model,  the  distribution  evolution  functions  Eqs.  (la)  and  (lb)  can 
lead  to  the  Navier-Stokes  equations  for  two-phase  system  [37]  and 
an  order  parameter  equation  for  phase  interface  evolution  under 
the  low  Mach  number  limitation: 

+  V  •  (pu)  =  0  (21a) 


+  V  •  (puu)  =  V  •  Pth  +  V  ■  (pi/Vu)  +  F  (21b) 

dt 

^+V-W  =  MV2/i  (21c) 

dt 

The  fluid  viscosity  v  is  related  to  the  relaxation  parameters  zy  in  the 
LBM  model  by 


[tf  -  0.5)5t 

v  =  — — ^ -  (22) 

In  the  unequal  viscosities  binary  liquid  system,  the  viscosity  at 
the  phase  interface  can  be  evaluated  by 


v  = 


(p  +  (po 
2ipo 


(v\  -Vg)  +  Vg 


(23) 


where  V\  and  vg  denote  the  viscosities  of  liquid  and  gas  phases 
with  the  equilibrium  order  parameter  of  ip0  and  -ip0,  respectively. 
The  density  p  is  reflected  in  the  order  parameter  ip  and  the  real 
local  density  can  be  obtained  in  the  terms  of  order  parameter  by 

P  =  Pg  +  (.Pi~  Pg)(V  +  lPo)R‘Po  [37], 

It  should  be  noted  that  the  surface  energy  between  the  solid 
wall  and  fluid  is  described  by  the  second  term  in  Eq.  (8)  with  S 
being  the  wall  surface  area.  Following  the  analysis  in  [42],  if  VK^s) 
is  expanded  as  a  power  series  in  the  wall  order  parameter  <ps  and 
keeping  only  the  first-order  term,  it  can  be  written  that  ^(<ps)  =  Xips 


The  velocity  or  pressure  boundary  condition,  which  is  related 
to  the  distribution  function /a(x,  t)  in  LBM,  can  be  dealt  with  using 
the  bounce-back  of  the  non-equilibrium  distribution  rule  as  given 
by  Hao  and  Cheng  [44]  for  the  D3Q19  model.  Moreover,  the  order 
parameter  boundary  for  the  distribution  function  g«(x,  t)  at  inlet  or 
outlet  is  an  additional  boundary  condition  in  the  two-phase  LBM 
compared  with  the  single-phase  model.  In  this  study,  the  pure  gas 
fluid  component  boundary  condition  is  used  at  both  the  inlet  and 
outlet  of  the  gas  channel.  This  boundary  condition  can  be  dealt  with 
similar  to  the  implementation  of  concentration  boundary  condi¬ 
tions  for  diffusion  problem  [45].  For  example,  if  the  inlet  boundary 
face  is  assumed  to  be  perpendicular  to  the  x-direction  with  the  lat¬ 
tice  velocity  e\ ,  e7,  eio,  e\\  and  e u  (as  shown  in  Fig.  2)  pointing  into 
the  calculation  region,  the  distribution  functions  g\ ,  g7,  gi0,  gu  and 
gi4  are  not  determined  after  streaming.  In  order  to  ensure  the  order 
parameter  (pin  at  inlet  equal  to  the  prescribed  value,  these  unknown 
distribution  functions  must  satisfy, 

gl  +  §7  +  glO  +  Su  +  &14  =  (pin  ~  (#0  +  g2  +  g3  +  g4  +  g5  +  g6  +  g8 


+  g9  +gl2  +gl3  +gl5  +gl6  +gl7  +gl8) 


(28) 


according  to  Eq.  (5).  Assuming  that  gt  (i  =  1,  7,  10,  11,  14)  are  dis¬ 
tributed  by  their  weight  coefficients  wa  given  in  Eq.  (3),  it  results 
in, 


Wiip* 

Wi  +  W7  +  Wio  +  Wn  +  Wi4 


i  =  1,  7,  10,  11,  14. 


(29) 


where  ip*  is  the  right-hand  side  of  Eq.  (28). 

In  order  to  implement  the  wetting  boundary  condition  at  the 
solid  wall  of  a  two-phase  system,  the  method  proposed  by  Niu  et 
al.  [37]  is  adopted  in  the  present  study.  In  this  method,  the  order 
parameter  derivative  in  Eq.  (24)  is  calculated  by  the  first-difference 
as  dip/dh  =  ( tpf  -  ips)8x  with  <ps  being  the  order  parameter  of  the 
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solid  and  cpf  the  order  parameter  of  the  fluid  point  adjacent  to  the 
solid  wall.  By  substituting  the  differences  into  Eq.  (24)  and  aver¬ 
aging  them  over  the  all  fluid  points  adjacent  to  the  solid  wall,  the 
order  parameter  cps  can  be  approximated  by 


<Ps  = 


N 


(30) 


where  N  is  the  total  number  of  the  fluid  points  which  are  nearest  the 
solid  wall.  Eq.  (30)  can  easily  be  applied  to  complex  solid  boundary 
such  as  porous  media.  As  this  is  our  first  numerical  simulation  of 
droplet  dynamics  in  gas  flow  channel  using  the  free-energy  method, 
we  limited  the  GDL  surface  to  be  smooth  in  present  study.  This 
model  can  easily  be  extended  to  study  the  roughness  effect  of  the 
fibrous  GDL  surface  on  droplet  dynamics. 

In  the  free-energy  LBM  model,  the  computation  of  derivatives 
vcp  and  v2(p  are  needed  to  evaluate  the  chemical  potential  given  by 
Eq.  ( 10)  and  pressure  tensor  given  by  Eq.  ( 14)  as  well  as  equilibrium 
distribution  functions  given  by  Eqs.  (2a)  and  (2b).  To  minimize  the 
discretization  error,  these  derivatives  are  calculated  using  19-point 
finite-difference  stencils  as  follows  [32,37]: 


oW“^r  +  e«St)e<*k 


Vk<p(r)  = 


StEa=0W«elk 


V2<p(r )  = 


2Z)«=ow“[<°(r  +  eaSt)  -  y(r)] 
Ealo  w«(eakSt)2 
with  k{k=x,y,  z )  denoting  the  dimension. 


(31a) 

(31b) 


Fig.  3.  Comparison  of  the  LBM  results  with  the  Laplace  law  for  pressure  jump  across 
a  stationary  droplet  interface. 


range  of  contact  angle  from  40°  to  160°  including  the  results  pre¬ 
sented  in  Fig.  4.  Therefore,  the  present  LBM  model  is  capable  of 
simulating  two-phase  flow  with  the  presence  of  solid  walls  accu¬ 
rately. 

5.  Results  and  discussion 


4.  Validation  of  the  present  model 

To  validate  the  present  LBM  model,  two  examples  of  a  droplet 
placing  in  an  unbounded  gas  domain  and  on  a  solid  wall  are  simu¬ 
lated.  For  the  first  benchmark,  a  two-dimensional  cylindrical  liquid 
droplet  is  initially  located  at  the  center  of  the  lattice  domain  with 
100  x  100  lattices  in  the  xy-plane.  The  periodic  boundary  condi¬ 
tions  are  imposed  at  all  boundaries.  According  to  the  Laplace  law, 
when  the  system  reaches  the  equilibrium  condition,  the  pressure 
difference  between  the  interior  and  exterior  of  the  droplet  A p  is 
related  to  the  surface  tension  a  as 

A  p=j  (32) 

where  R  is  the  droplet  radius.  Fig.  3  shows  the  comparison  of  simu¬ 
lated  pressure  jump  across  the  droplet  interface  versus  the  inverse 
of  droplet  radii  with  the  Laplace  law  given  by  Eq.  (32)  for  a  given 
surface  tension  of  0.2.  As  shown,  the  simulation  results  are  in  excel¬ 
lent  agreement  with  the  Laplace  law,  which  is  represented  by  a  solid 
straight  line. 

Next,  the  equilibrium  contact  angle  of  a  droplet  placed  on  the 
solid  wall  is  simulated  in  a  computational  domain  of  60  x  60  x  50 
lattices.  The  periodic  boundaries  are  used  in  x-  and  y-directions 
while  wall  boundaries  at  top  and  bottom  in  z-direction  are  imposed. 
Initially,  a  water  droplet  of  a  cubic  shape  is  located  on  the  bot¬ 
tom  wall,  and  the  simulation  is  performed  until  the  droplet  reaches 
its  equilibrium  with  an  unchanged  spherical-cap  shape.  The  cor¬ 
responding  equilibrium  contact  angle  of  the  droplet  on  the  wall 
is  obtained  by  prescribing  a  certain  wetting  potential  0  according 
to  Eq.  (26).  Fig.  4(a)-(c)  shows  equilibrium  shapes  of  the  droplet 
with  the  wetting  potential  0  being  0.3,  0  and  -0.3,  respectively. 
Their  corresponding  equilibrium  contact  angles,  obtained  from 
the  interface  topologies  are  61°,  89°  and  118°,  respectively.  The 
simulated  equilibrium  contact  angle  as  a  function  of  the  wet¬ 
ting  potential  for  a  wider  range  is  presented  in  Fig.  5,  where  it 
is  shown  that  the  simulated  results  by  the  present  LBM  model 
agree  quite  well  with  the  theoretical  predictions  by  Eq.  (26)  in  the 


5.1.  Droplet  formation  in  gas  flow  channel 

Fig.  1(b)  shows  the  schematic  of  a  droplet  emergence  through 
the  GDL.  The  computational  domain  includes  a  micro-channel  with 
the  cross  section  of  60  (width)  x  30  (height)  lattices  and  the  length 
of  120  lattices,  as  well  as  a  GDL  of  10  lattices  thickness  with  a  sin¬ 
gle  square  pore  locating  at  the  midline  of  the  channel.  The  side 
length  of  the  micro-pore  is  9  lattices  and  one  lattice  distance  rep¬ 
resents  10  |jim  of  the  physical  length  in  present  study.  Thus,  the 
micro-channel  has  a  representative  geometrical  size  of  gas  flow 
channel  used  in  a  micro-PEMFC  [46]  and  the  micro-emergence 
pore  has  its  diameter  close  to  the  bigger  pores  of  the  carbon  paper 
GDL  [20].  This  is  consistent  with  experimental  observation  that 
water  usually  emerges  on  the  GDL  surface  along  the  fixed  paths 
including  larger  and  less  hydrophobic  pores  in  GDL  [8,17-19,47]. 
The  contact  angle  of  the  pure  PTFE  material  is  reported  to  be 
~110°  and  the  contact  angle  of  the  GDL  surface  treated  by  PTFE 
can  reach  to  170°  due  to  roughness  [19].  Therefore,  the  contact 
angle  of  the  bottom  hydrophobic  GDL  surface  is  considered  to 
be  in  the  range  of  110-170°  in  this  study.  Additionally,  since  the 
interaction  between  the  droplet  and  the  sidewalls  is  not  consid¬ 
ered  here,  it  is  assumed  that  the  wettability  of  the  two  sidewalls 
and  top  wall  are  all  hydrophilic  with  60°  contact  angle.  In  addi¬ 
tion,  the  velocity  of  the  water  through  the  square  pore  is  fixed  at 
0.075 ms-1,  which  is  close  to  the  value  used  by  Chen  et  al.  [15]. 
The  Reynolds  number  Re  for  gas  flow  is  restricted  to  be  smaller 
than  100  with  the  channel  equivalent  diameter  (400  pan)  being 
the  characteristic  length.  To  relate  the  physical  space  with  lattice 
space,  a  length  scale  l0  =  10|xm  and  a  time  scale  t0  =  3.5  x  10-7  s 
are  chosen,  which  results  in  the  physical  velocity  up  rescaled  to 
the  lattice  velocity  14  by  i4  =  upto/fo-  The  grid  effect  is  tested  by 
performing  additional  simulations  in  a  finer  grid  with  /0  =  5  pm. 
The  results  show  similar  behavior  of  water  droplet  dynamics  for 
both  two  grid  sizes  and  the  difference  of  the  formed  droplet  size 
between  them  is  below  5%.  So  the  grid  with  /0  =  10  [xm  is  acceptable 
in  this  study.  The  simulated  results  by  the  LBM  model  on  dynamics 
of  the  droplet  emerging  and  moving  in  the  gas  flow  channel  with 
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Fig.  4.  Droplets  on  a  surface  with  the  wetting  potential  (9  of  0.3  (a),  0  (b)  and  -0.3 
(c).  The  corresponding  contact  angles  are  61°,  90°  and  118°,  respectively. 


various  gas  entering  velocities  and  GDL  wettability  are  presented 
below. 

Fig.  6  shows  the  snapshots  of  droplet  formation  simulations  for 
the  bottom  wall  equilibrium  contact  angle  of  120°  and  gas  side  cap¬ 
illary  number  Ca  (defined  as  Ca  =  riaUinla,  where  Uin  is  the  average 
inlet  gas  velocity,  r\a  is  the  gas  viscosity  and  o  is  the  surface  tension) 
of  0.06.  Due  to  the  dominate  effect  of  surface  tension,  the  emerging 


Fig.  5.  Comparison  of  the  simulated  equilibrium  contact  angles  with  the  theoretical 
predictions  by  Eq.  (26). 


water  droplet  keeps  an  approximate  spherical-cap  shape  above  the 
emergence  pore  (see  Fig.  6(a)).  As  the  growth  of  the  droplet,  the  drag 
force  exerted  on  the  droplet  from  gas  shear  increases  and  forces  the 
droplet  to  creep  downstream  slowly  (Fig.  6(b)).  Flowever,  before  the 
drag  force  overcomes  the  resistance  which  is  dominated  by  the  sur¬ 
face  tension  force  due  to  the  connection  between  the  droplet  and 
the  micro-pore  at  this  stage,  the  liquid  water  connecting  bridge 
is  maintained  with  the  incessant  emergence  of  water  through  the 
pore,  and  is  elongated  gradually  at  the  same  time  (Fig.  6(c)),  result¬ 
ing  in  the  formation  of  a  neck  after  the  droplet  (Fig.  6(d)).  Once  the 
droplet  grows  larger  and  when  the  shearing  stress  on  the  droplet 
exceeds  the  surface  tension  force,  it  will  accelerate  and  move  down¬ 
stream  and  thus  the  neck  ruptures,  inducing  the  final  detachment 
of  the  droplet  (Fig.  6(e)).  Subsequently,  the  detached  droplet  moves 
to  the  outlet  of  the  channel  along  the  bottom  wall  and  simultane¬ 
ously  a  fresh  droplet  is  formed  through  the  micro-pore  repeatedly 
as  shown  in  Fig.  6(f).  The  simulation  results  in  Fig.  6  capture  the 
dynamic  evolution  of  the  droplet  formation  and  its  subsequent 
motion  on  the  GDL  surface,  which  are  accordant  with  experimental 
visualizations  [14].  From  the  above  analysis,  it  can  be  concluded  that 
the  mechanism  of  the  water  droplet  detachment  on  the  assumed 
GDL  surface  can  be  regarded  as  the  competition  between  the  sur¬ 
face  tension  force  from  the  connecting  of  the  water  in  the  GDL  pore 
and  the  emerged  part  on  the  GDL  surface,  and  the  drag  force  which 
is  related  to  the  gas  flow  velocity  and  the  droplet  size  in  the  chan¬ 
nel.  Similar  conclusions  were  also  been  reported  by  Kimball  et  al. 
in  their  most  recent  experimental  studies  [17,18].  Fig.  7  shows  the 
enlarged  view  of  the  profile  of  the  detached  droplet  in  Fig.  6(f)  and 
the  velocity  field  relative  to  the  droplet.  The  deformation  of  the 
moving  droplet  and  its  rolling  motion  on  the  hydrophobic  surface 
is  clearly  shown  in  this  figure. 

The  effect  of  the  GDL  surface  wettability  on  droplet  dynamics 
can  be  shown  by  performing  simulations  with  different  equilib¬ 
rium  contact  angles.  Fig.  8  shows  the  droplet  evolution  process  with 
bottom  wall  equilibrium  contact  angle  of  165°  and  the  capillary 
number  of  0.06.  Compared  with  Fig.  6  with  wall  equilibrium  con¬ 
tact  angle  of  120°,  the  more  hydrophobic  surface  causes  a  droplet 
with  the  larger  height  with  the  same  water  volume  (Fig.  8(a) 
and  (b)),  thus  leading  to  a  larger  drag  force  on  it.  As  a  result, 
the  appearance  of  the  neck  (Fig.  8(c))  and  the  detachment  of  the 
droplet  (Fig.  8(d))  occur  earlier.  Flowever,  after  the  detachment,  the 
droplet  seems  to  be  lifted  from  the  bottom  surface  with  the  con¬ 
tact  area  decreasing  gradually  as  it  moves  towards  to  the  channel 
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(a)  t=  1.12  ms 


(d)/  =  4.76  ms 


(b)  t  =  2.80  ms 


(e)  t  =  5.18  ms 


(c)f  =  4.20ms  (f)/  =  6.65ms 


Fig.  6.  Snapshots  of  the  droplet  formation  (blue  area  denotes  the  droplet)  through  a  micro-pore  in  the  bottom  wall  with  the  equilibrium  contact  angle  of  120°  and  Ca  =  0.06. 
The  upper  2D  figures  are  the  cross-sectional  view  on  the  middle  section  along  the  gas  flow  channel.  (For  interpretation  of  the  references  to  color  in  the  figure  caption,  the 
reader  is  referred  to  the  web  version  of  the  article.) 


outlet  (Fig.  8(e)).  This  can  be  attributed  to  the  larger  lift  force  than 
adhesion  force  acting  on  the  spherical  droplet  in  the  channel.  In 
general,  the  adhesion  force  of  partially  wetting  liquid  on  the  solid 
wall  is  due  to  surface  tension  and  this  force  can  be  approximated 
expressed  by  Fadhension  =  2ttRg  sin2  0  for  a  spherical-cap  droplet.  If 
the  gravity  is  neglected  for  a  micro-droplet  here,  the  droplet  would 


lift  from  the  surface  when  the  lift  force  overcomes  this  adhesion 
force  (Fig.  8(f)).  The  lift  force  is  a  function  of  the  gas  flow  veloc¬ 
ity  [48,49]  and  increases  with  the  gas  flow  velocity.  Therefore, 
it  can  be  concluded  that  the  more  hydrophobic  surface  and  the 
larger  gas  flow  velocity  mean  smaller  adhesion  force  and  bigger 
lift  force,  thus  are  beneficial  for  the  droplet  to  lift  from  the  GDL  sur- 
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Fig.  7.  Enlarged  view  of  the  profile  of  the  moving  droplet  in  Fig.  6(f)  and  the  corre¬ 
sponding  velocity  field.  The  velocity  field  is  relative  to  the  droplet  velocity. 

face.  Consequently,  more  GDL  surface  will  be  available  for  reactant 
transport. 

The  dependence  of  the  detached  droplet  diameter  ( Ddwpiet ) 
on  the  gas  capillary  number  Ca  and  equilibrium  contact  angle 
0  is  shown  in  Fig.  9(a)  and  (b),  respectively.  As  expected,  the 
droplet  detachment  size  decreases  with  the  increase  of  Ca  and  GDL 
hydrophobicity.  However,  their  influence  becomes  weaker  when  Ca 
or  0  is  higher  than  a  certain  value.  The  simulation  results  of  droplet 
dynamic  behavior  in  gas  flow  channel  in  present  study  are  simi¬ 
lar  to  the  results  obtained  in  most  recent  VOF  studies  [23-26].  As 
mentioned  above,  the  LBM  model  has  the  advantages  of  no  need  to 
track  phase  interface  and  the  mesoscopic  insight  of  contact  angle. 
Thus,  it  is  a  promising  method  to  investigate  multiphase  flow  such 
as  water  transport  in  PEMFCs. 


the  highly  hydrophobic  surface,  resulting  in  the  little  influence  of 
droplet  size  on  its  moving  velocity. 

5.3.  Analytical  model  on  emerging  droplet  in  micro-channel 

Following  the  work  by  Chen  et  al.  [15]  as  well  as  by  Kumbur 
et  al.  [16],  we  now  carry  out  a  similar  analysis  on  an  emerging 
droplet  from  a  micro-pore  according  to  force  balance  considera¬ 
tion.  We  consider  the  specific  case  of  a  micro-gas  flow  channel  with 
a  hydrophobic  GDL  surface  because  of  the  usual  PTFE  treatment. 
Since  surface  tension  plays  a  dominant  effect  on  droplet  behaviors 
in  micro-channels,  the  capillary  number  is  a  particularly  important 
dimensionless  parameter  in  this  study.  Moreover,  the  Bond  number 
Bo  (defined  as  Bo  =  pgD2lcr  with  g  being  the  gravity  acceleration) 
is  much  smaller  than  unity  in  present  problem,  implying  that  the 
gravitational  force  can  be  neglected.  Furthermore,  in  order  to  sim¬ 
plify  the  analysis,  other  assumptions  are  made  as  follows:  (i)  the  gas 
flow  in  the  micro-channel  is  fully  developed  laminar  flow;  (ii)  the 
width  of  the  channel  is  larger  than  its  height,  so  that  the  sides  of  the 
channel  have  small  influence  on  the  gas  flow;  (iii)  the  droplet  keeps 
a  spherical-cap  approximately  with  the  neglecting  of  the  devia¬ 
tion  from  the  spherical  shape  because  of  the  small  viscous  force.  It 
should  be  noted  that  similar  force  analysis  models  have  been  used 
to  investigate  the  liquid  water  droplet  instability  on  a  GDL  surface 
[15,16,50].  In  these  studies,  the  droplet  instability  is  controlled  by 
the  surface  tension  force  induced  by  the  contact  angle  hysteresis. 
However,  in  order  to  correspond  to  the  above  simulation  condition 
for  a  smooth  and  homogenous  surface,  the  static  contact  angle  hys¬ 
teresis  of  a  droplet  is  neglected  in  our  analytical  model.  Thus,  the 
dominant  surface  tension  force  in  this  study  is  different  from  that 
in  droplet  instability  studies. 

As  shown  in  Fig.  1(b),  the  control  volume  is  denoted  by  the 
dashed  line  enclosing  the  whole  droplet  and  the  volume  above  it 
[15,16].  We  choose  the  frame  of  reference  with  x,  y-axes  parallel 
and  perpendicular  to  the  wall,  respectively.  The  origin  of  y-axis  is 
fixed  on  the  top  wall.  According  to  the  spherical  shape  assumption, 
the  droplet  height  h  and  the  droplet  contact  diameter  l  on  the  GDL 
surface,  as  well  as  the  gas  flow  width  b  above  the  droplet  can  be 
obtained  from  the  droplet  and  channel  geometry  as  follows: 


5.2.  Droplet  removal  time 


h  =  f?(l  -  cos  0)\  l  =  2ft sin 0;  b  =  H-h  (33) 


The  time  of  the  water  droplet  removing  out  of  the  gas  flow  chan¬ 
nel  is  a  key  parameter  in  PEMFCs  water  management.  It  is  expected 
that  the  total  time  of  water  removal  is  consisting  of  two  stages, 
namely  the  formation  stage  and  the  moving  stage.  The  formation 
time  tf  can  approximately  be  obtained  from  the  flow  rate  Qd  of 
the  emerging  liquid  water  and  the  size  of  the  detached  droplet 
by  tf=VdIQd  with  Vd  being  the  detached  droplet  volume,  while 
the  moving  time  tm  is  expected  to  be  proportional  to  the  mov¬ 
ing  channel  length,  which  is  considered  to  be  the  distance  from 
the  emergence  pore  to  the  channel  outlet  (100  |jim).  The  simulated 
formation  time  tf  and  total  removal  time  ( tf+  tm)  of  the  droplet  at 
different  capillary  numbers  Ca  and  GDL  equilibrium  contact  angles 
are  presented  in  Fig.  10(a)  and  (b),  respectively.  As  shown  from 
the  figures,  larger  Ca  and  GDL  equilibrium  contact  angle  result  in 
the  quicker  removal  of  the  water  droplet  away.  However,  the  total 
removal  time  becomes  close  to  the  formation  time  gradually  with 
the  increase  of  the  gas  flow  velocity  (i.e.,  larger  value  of  Ca),  which 
implies  that  the  water  droplet  emergence  process  will  dominate 
the  water  droplet  removal  time  under  large  gas  flow  velocity  con¬ 
ditions.  Moreover,  the  benefit  of  increasing  hydrophobicity  of  GDL 
larger  than  145°  becomes  unobvious  for  water  droplet  removal,  as 
the  total  removal  time  keeps  fairly  flat  when  the  equilibrium  con¬ 
tact  angle  beyond  145°,  although  the  size  of  the  formed  droplet  is 
still  decreasing.  This  is  probably  due  to  the  lifting  of  the  droplet  from 


where  H  is  the  channel  height,  R  is  the  droplet  radius  and  0  is  the 
equilibrium  contact  angle.  For  a  fully  developed  laminar  flow,  the 
gas  flow  profile  above  the  droplet  can  be  assumed  to  be  of  the 
Poiseuille  type,  which  satisfies  the  Stokes’  equation: 


d2u  _  dp 
=  dx 


(34) 


where  pa  denotes  the  gas  viscosity.  In  addition,  the  moving  velocity 
of  the  emerging  droplet  is  ignored  here.  This  assumption  is  approx¬ 
imately  satisfied  due  to  the  low  creeping  velocity  of  the  droplet  at 
its  growth  stage  and  has  been  proved  by  simulations  above.  As  a 
result,  the  solution  of  Eq.  (34)  can  be  obtained  as 


6Ua  2 

u=~wy  + 


dp  _  6  Ua 

-  -2ria~tiT 


(35) 


where  Ua  is  the  average  gas  velocity  along  the  flow  direction  above 
the  droplet  and  is  related  to  the  average  inlet  gas  velocity  Uin  by 


Ua  =  "uir 


(36) 


For  a  Newtonian  fluid,  the  gas  shear  stress  at  the  top  wall  can  be 
calculated  by 


c  _  9ui  c 

Tsftear  —  py  '^a  ~  ha 


6 Ua 

b 


(37) 
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(b)  /  =  2.80  ms  (e)  t  =  4.90  ms 


(c)  t  =  3.50  ms  (f)  t  =  6. 1 6  ms 

Fig.  8.  Snapshots  of  the  droplet  formation  (blue  area  denotes  the  droplet)  through  a  micro-pore  in  the  bottom  wall  with  the  equilibrium  contact  angle  of  165°  and  Ca  =  0.06. 
The  upper  2D  figures  are  the  cross-sectional  view  on  the  middle  section  along  the  gas  flow  channel.  (For  interpretation  of  the  references  to  color  in  the  figure  caption,  the 
reader  is  referred  to  the  web  version  of  the  article.) 


where  Sa  =  {2R )2  is  the  contact  area  of  the  control  volume  on  the 
top  wall.  The  pressure  force  Fp  along  the  flow  direction  is  related 
to  the  pressure  gradient  given  by  Eq.  (35)  and  the  distance  over  the 
droplet  by 


with  Sp  =  2 HR  being  the  cross  sectional  area,  which  is  perpendicular 
to  the  flow  direction.  Substituting  Eq.  (36)  into  Eqs.  (37)  and  (38), 
the  drag  force  acting  on  the  droplet  in  the  gas  flow  channel  can  be 
obtained  from  the  combination  of  the  two  terms  by 


Fp  —  ( Pi  —  Po)Sp  —  ^Pa-j^-  •  2 RSp 


(38) 


_  TT  48 H2R2  tt  24 HR2 

Fdrag  ~  PaClin  ^3  Pamiri 


(39) 
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Fig.  9.  (a)  Comparison  of  the  LBM  simulated  droplet  diameters  with  the  predictions 
by  Eq.  (42)  for  different  capillary  numbers  Ca  at  equilibrium  contact  angle  of  120°. 
(b)  Comparison  of  the  LBM  simulated  droplet  diameters  with  the  predictions  by  Eq. 
(42)  for  different  equilibrium  contact  angles  6  at  Ca  =  0.06. 


100  110  120  130  140  150  160  170 

Equilibrium  Contact  Angle  (°) 

Fig.  10.  (a)  The  simulated  formation  and  removal  time  of  the  droplet  for  different 
capillary  numbers  Ca  at  equilibrium  contact  angle  of  120°.  (b)  The  simulated  forma¬ 
tion  and  removal  time  of  the  droplet  for  different  equilibrium  contact  angles  6  at 
Ca  =  0.06. 


Besides,  the  resistance  force  Fsurface  is  directly  related  to  the  sur¬ 
face  tension  and  can  further  be  divided  into  two  terms:  the  surface 
tension  force  due  to  the  connecting  of  the  droplet  to  the  emergence 
pore,  and  the  surface  tension  force  due  to  the  droplet  deforma¬ 
tion.  The  first  term  is  related  to  the  size  of  the  emergence  pore 
as 

Fsurface,  1  =  ~Cf7td  (40a) 

with  d  being  the  equivalent  diameter  of  the  emergence  pore  and 
o  the  surface  tension  of  the  liquid/gas  interface.  It  should  be  noted 
that  this  surface  tension  force  term  was  not  considered  in  previ¬ 
ous  investigations  on  stable  characteristic  of  a  static  droplet  [15,16], 
although  it  has  significant  influence  on  the  droplet  dynamics  and 
formation  [20]  and  has  been  proved  to  be  the  dominant  cohesive 
force  for  the  detachment  of  the  emerged  droplet  by  experiments 
[18].  The  second  term  of  resistance  force  is  related  to  the  difference 
between  the  advancing  and  receding  contact  angles.  If  the  contact 
line  on  the  wall  is  assumed  to  be  a  circle  and  its  semicircle  at  the 
front  or  rear  part  of  the  contact  line  is  assumed  to  be  constant 
advancing  or  receding  contact  angle,  respectively,  this  resistance 
force  is  given  by  [51]: 

Fsurface, 2  =  -<7(C0S  6r  -  COS  9a)L  (40b) 


where  6r  and  0a  are  the  receding  contact  angle  and  the  advancing 
contact  angle  as  shown  in  Fig.  1(b),  and  L  is  the  contact  line  length 
of  the  front  or  rear  semicircle  part,  given  by  L  =  tt// 2.  For  a  rough 
surface  with  a  non-moving  droplet,  the  term  Eq.  (40b)  is  related  to 
the  contact  angle  hysteresis,  which  is  a  function  of  surface  proper¬ 
ties  such  as  roughness.  However,  for  a  smooth  surface  with  contact 
line  motion,  it  is  related  to  the  dynamic  contact  angle,  which  has 
been  shown  to  be  a  function  of  the  dynamic  contact  line  velocity 
[52-55]. 

Finally,  the  force  balance  along  the  x-direction  on  the  control 
volume,  which  includes  an  emerging  droplet,  results  in: 


Fp  +  Fshear  +  F surface  ~  0  (41 ) 

Therefore,  Eq.  (41 )  provides  a  critical  criterion  for  the  formed 
droplet  detaching  from  the  emerging  micro-pore.  Additionally,  due 
to  the  much  smaller  value  of  the  droplet  creeping  velocity  than  the 
gas  flow  velocity,  the  force  Fsurface2  deriving  from  droplet  motion 
is  regarded  to  be  much  smaller  than  Fsurface tl ,  thus  can  be  ignored 
here.  As  a  result,  substituting  Eqs.  (37),  (38)  and  (40a)  into  Eq.  (41 ), 
the  force  balance  of  an  emerging  droplet  can  further  be  reduced  to: 


48 


H2R2 

-fi¬ 


nd 

Ca 


(42) 
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So,  if  Eq.  (33)  is  substituted  into  Eq.  (42),  the  size  of  the  detached 
droplet  can  be  predicted  in  a  given  channel.  The  corresponding  pre¬ 
diction  results  of  the  droplet  size  as  a  function  of  Ca  and  equilibrium 
contact  angle  0  are  also  plotted  in  Fig.  9(a)  and  (b)  by  a  solid  line, 
respectively.  As  shown,  the  predicted  values  are  slightly  smaller 
than  the  simulation  results.  This  underestimation  of  droplet  size  can 
be  attributed  to  the  neglect  of  droplet  creeping  velocity  in  the  ana¬ 
lytical  model,  resulting  in  the  underestimation  of  resistance  force. 
However,  theoretical  prediction  of  this  force  given  by  Eq.  (40b) 
has  some  difficulties  because  of  the  unclear  physical  mechanism 
of  dynamic  contact  angle  with  contact  line  motion.  Fortunately,  the 
LBM  model  is  capable  of  simulating  the  change  of  the  dynamic  con¬ 
tact  angle  at  mesoscopic  level,  thus  provides  more  accurate  results 
in  this  study. 

6.  Concluding  remarks 

In  this  work,  the  dynamic  behavior  of  water  droplet  formation 
and  removal  in  the  micro-gas  flow  channel  of  a  PEMFC  has  been 
simulated  using  the  multiphase  free-energy  LBM  approach.  The 
evolution  process  of  water  emerging  through  a  micro-pore  is  clearly 
captured  by  simulation  and  is  shown  in  good  agreement  with  visu¬ 
alization  results.  The  size  of  the  detached  droplet  and  the  time  of 
the  liquid  water  removing  out  of  the  channel  are  investigated  with 
various  gas  flow  velocity  and  GDL  surface  wettability.  The  results 
show  that  the  increase  of  the  gas  flow  velocity  and  the  enhancement 
of  the  GDL  hydrophobicity  can  facilitate  the  formation  of  smaller 
droplets,  thus  decreasing  the  time  of  water  droplets  adhering  to 
the  GDL  surface  near  the  emergence  pore.  In  addition,  the  highly 
hydrophobic  surface  is  shown  to  promote  droplet  lifting  from  the 
surface  after  its  detachment,  resulting  in  larger  GDL  surface  avail¬ 
able  for  gas  reactant  transport.  An  analytical  model  is  also  presented 
to  predict  the  droplet  departure  sizes  based  on  a  force  balance  con¬ 
sideration,  and  results  from  this  analysis  are  in  good  agreement 
with  that  from  LBM  simulation  under  various  conditions. 

The  results  of  the  present  study  provide  a  mesoscopic  insight  for 
understanding  the  water  droplet  dynamic  behavior  in  the  gas  flow 
channel  of  PEMFCs,  which  is  helpful  for  solving  the  channel  flood¬ 
ing  problem.  However,  it  should  be  noted  that  the  static  contact 
angle  hysteresis,  a  typical  characteristic  of  an  actual  roughness  sur¬ 
face  such  as  GDL  surface,  is  not  considered  in  this  work.  Therefore, 
further  studies  on  this  aspect  of  the  problem  are  needed. 
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